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1 Introduction 



Quantum chromodynamics (QCD) has been very successful in describing the high- 
energy interactions of quarks and gluons; however, despite nearly two decades of 
lattice simulations, attempts to extract from the theory reliable quantitative infor- 
mation on the low-energy properties of hadrons have not been fully successful. 

Recently, a new approach has been suggested M which advocates combining effec- 
tive field theory and lattice techniques to study hadrons composed entirely of heavy 
quarks. There are many good reasons for doing this. Such hadrons are the simplest to 
analyze since the quarks are nonrelativistic. Simulations which use a nonrelativistic 
QCD (NRQCD) action instead of the usual Dirac action are far more efficient since 
the heavy-quark propagator can be computed as an initial-value problem in a single 
sweep through the lattice. The Dirac propagator is a boundary- value problem which 
must be solved iteratively using many sweeps. Also, the fermion doubling problem 
does not occur in the nonrelativistic theory. Heavy quarkonia, such as the ip and T 
mesons, are small so that lattice volumes presently used are sufficiently large, and 
their properties are fairly insensitive to both light-quark vacuum polarization effects 
0, which are very expensive to simulate, and to heavy-quark vacuum polarization. 
Simulations of heavy quarkonia do not suffer from a rapid reduction in the signal-to- 
noise ratio which is a serious shortcoming of the static approximation used to study 
heavy-light systems, such as the D and B mesons. Quarkonia are well understood; 
this makes possible good control of systematic errors which is necessary for precise 
tests of QCD. Furthermore, an abundance of experimental data on quarkonium states 
is available to which the simulation results, such as level splittings, decay constants, 
and wavef unctions, may be compared. 

NRQCD is an effective cutoff field theory constructed from a set of nonrenormal- 
izable interactions specified solely by the symmetries of QCD, the chosen regulator, 
locality, and the accuracy desired. It is essentially a low-energy expansion of the 
Dirac theory in terms of the expectation value v of the heavy quark velocity in a typ- 
ical heavy quark hadron. An action which includes all spin-independent relativistic 
interactions suppressed by v 2 relative to the leading terms and all spin- dependent cor- 
rections up to order v 4 has previously been formulated using a lattice regularization 
scheme @ and is here referred to as lattice NRQCD. To fully define lattice NRQCD, 
the coupling coefficients of the interactions appearing in the action must be specified. 
These are process and momentum independent and are uniquely determined (for a 
given regulator) by requiring that lattice NRQCD exactly reproduces the results of 
continuum QCD at low energies. 

Since the role of these couplings is to absorb the relativistic effects arising from 
highly-ultraviolet QCD processes, one expects that they may be computed to a good 
approximation using perturbation theory, provided the quark mass M is large enough. 
The simplest way to proceed is to evaluate various scattering amplitudes both in 
QCD and lattice NRQCD and adjust the couplings until these amplitudes agree to 
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the desired order in v and the QCD coupling g. In this way, one obtains coupling 
coefficients which are power series in g 2 (A), where A is the cutoff of the effective 
theory. The cutoff A must be large enough so that 0(g 2 (A)) corrections to the 
effective couplings are small, making a perturbative analysis meaningful. However, 
power-law divergences generally occur, producing terms such as g 2 (A)A/M which 
render perturbation theory useless if A is made too large. 

In this paper, the lowest-order corrections to the heavy quark self-energy in lattice 
NRQCD are calculated using weak-coupling perturbation theory. The action formu- 
lated in Ref. p| which includes spin effects is used. Similar to an earlier study [Q, 
the mass and wavefunction renormalization parameters required to match continuum 
QCD are obtained, as well as a necessary overall energy shift. Lattice NRQCD is 
briefly described in Section The Feynman rules are derived in Section |3] and the 
self-energy calculations are presented and discussed in Section £|. A tadpole improve- 
ment scheme in which all link variables are rescaled by a mean-field factor (link vari- 
able renormalization) is also studied. Taking into account the mean-field corrections 
introduced by this scheme, the coefficients of g 2 in the heavy quark renormalization 
parameters are found to be small. Section |5| offers conclusions. 

2 Lattice NRQCD 

Standard renormalization group techniques are used to formulate the NRQCD action. 
First, since the physics of heavy quark systems is dominated by momenta p ~ Mv <C 
M, where M is the heavy quark mass, an ultraviolet cutoff A ~ M is introduced. The 
Dirac action can then be replaced by a nonrelativistic Schrodinger action in which 
the quark and antiquark degrees of freedom essentially decouple. However, new local 
interactions must be added in order to compensate for the loss of the discarded 
relativistic states. Power counting rules @ are used to classify these interactions 
according to their estimated magnitude in a typical quarkonium state as measured 
in terms of v. Equipped with these rules, the interaction terms which contribute to 
quarkonia physics up to some given order in v may then be enumerated. The couplings 
associated with these interactions are determined by requiring agreement between 
the physical results of NRQCD and those of full QCD through the specified order in 
the quark velocity. When formulated on a lattice, NRQCD becomes a particularly 
powerful tool for studying heavy quark systems. 

In lattice simulations, Euclidean space correlation functions are computed. These 
may be obtained from the Minkowski theory by Wick rotating the contour integral 
of the Lagrangian over time from the real axis to the imaginary axis in a clockwise 
manner. Also, the path integration over the scalar potential must be similarly Wick 
rotated at each point in space-time. Then by defining x 4 = ix°, £4 = — ix$, A 4 = iA°, 
and A 4 = —iA , one can again work with real quantities. The spatial components 
of all four-vectors remain unchanged. The path integral weight exp(iS) becomes 
exp(— Se), where Se is called the Euclidean action. Wick rotation transforms the 
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Minkowski metric tensor g^ = diag(l, —1, —1, —1) into the negative of the identity 
matrix, so that in Euclidean space, raising or lowering any index introduces a sign 
change. 

In lattice NRQCD, the quark fields are defined on the sites of a four- dimensional 
hypercubic lattice with spacing a, while the gauge-field degrees of freedom reside on 
the links between the sites. With each link originating at a site x and terminating 
at a site x + ae M is associated a link variable U^x), which is a lattice version of the 
parallel transport matrix between sites and is an element of the Lie group associated 
with the gauge invariance of the theory. A lattice gluon field A^(x) = Y%=i ^{ x )T b 
may be defined in terms of the link variables using 



U^x) = exp 



iagA^x + -e^) 



where g is the coupling constant of the theory and (e M ) a = 5°. This definition is 
convenient since it is simple and satisfies UUx) = LL^x+ae^). The SU(3) generators 
T b are Hermitian, traceless, and satisfy Tr(T a T fe ) = \5 ab and [T a , T b ] = if abc T c , where 
jabc are real, fully antisymmetric structure constants. Of course, the field A^(x) 
is not identical to the gluon field of the continuum theory. Let G>(x) represent the 
gluon field defined in the continuum theory. Under a local gauge transformation 
S(x), G>(x) — > S(x)G /Ji (x)S^(x) — (i/g)S(x)d ll S^(x). The field A^{x) does not 
transform in this way; rather, it transforms in a very complicated manner in order 
that the link variable transforms under a local gauge transformation according to 
U^ix) — > S(x)U fl (x)S 1 '(x + ae M ). 

In Minkowski space, the chromoelectric and chromomagnetic fields are usually 
defined in terms of the field strength tensor in the following manner: E? M \(x) = 
F k0 (x) = Fofc(x), and B k (x) = — ^€ki m F lm (x), where ei 23 = 1. In Euclidean space, 
the magnetic field is unchanged. The electric field in Euclidean space is here defined 
by E (x) = F k4 (x) = — F 4fc (x), so that E h M ^ = —iE k . The Hermitian and traceless 
field strength tensor is best represented 0] by cloverleaf operators defined at the sites 
of the lattice: 



F^{x) = T^ix) - -TrF^x), (2) 



•M*) = ^(n^xj-n^x)), (3) 

^fiu(x) = j U a (x)U^(x + ae a )U- a (x + ae a + aep)U-p(x + aej3), (4) 

{(":/3) V 

with {(a, f3)}fj, u = {(a*j v), {y, — /x), (— /x, — v), (— v, /x)} for /x ^ v. This representation 
is chosen since it transforms as the (1, 0) © (0, 1) six-dimensional reducible represen- 
tation of the hypercubic group, similar to the continuum case. 

Covariant derivatives are replaced by forward, backward, or symmetric differences 
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on the lattice: 



5 A« +) (!/;x) = U,(y)S^\x,y + ae,)-S^(x,y), (5) 



a 5 AH(y;x) = 5^(x,y)-Ul(x)5^(x,y-ae ti ), (6) 



= i(AW+AH), (7) 

where denotes a four- dimensional Kronecker ^-function. For example, 

E^fc^W = U^yMy+aeJ-^y). (8) 
Also, the Laplacian becomes 

a (2) = Ea1 +) aH = E^ +) . (9) 

fe=i fe=i 

The finite lattice spacing introduces systematic errors into NRQCD which can 
be reduced by the addition of new interactions to the action. At tree level, this is 
most easily accomplished by improving the components comprising the lattice action 
so that they more accurately reproduce the effects of their continuum counterparts. 
An improved difference operator which reproduces the behavior of the continuum 
covariant derivative through order a 4 is given by 



= m 



An improved lattice Laplacian is 



12 k=1 



Lastly, an improved cloverleaf field strength tensor may be defined by 



5 „ , , 1 r 



3 / 6 



x 



+ U^(x - ae^F^x - ae^U^x - ae^) - (/i <-> u) . (12) 

The portion of the lattice NRQCD action containing the heavy quark-gluon in- 
teractions may be written 

4 n) = a3 E^(^W (1 3 ) 



-a 



This action includes all spin-independent corrections which are suppressed by v 2 rel- 
ative to the leading terms and all spin-dependent interactions suppressed by factors 
up to v 4 . The heavy quark field i/j(x) is a Pauli spinor corresponding to the two upper 
components of the original Dirac field. (An analogous action may be written for the 
heavy antiquark field x( x )-) The improved kinetic energy is given by 

A( 2 ) a (A^) 2 , . 

#o = jt-, 14 

2M An AM 2 ' v ' 

where n is a positive integer. This operator is cleverly constructed to produce a leading 

error from the lattice approximation of the temporal derivative which can be removed 

by a redefinition of the quark fields. In this way, the heavy quark propagation remains 

governed by a Schrodinger equation across a single time slice without affecting energy 

levels, on-shell scattering amplitudes and other physical quantities. The parameter n 

was introduced in Ref. M to remove instabilities in the evolution of the quark Green's 

function which occur when the temporal spacing is not small enough to accurately 

treat the high-momenta modes. 

The quark-gluon interactions are as follows: 

5H = j2c 1 V j , (15) 



where 



v ^ = -W-' (l6) 

V 2 = ^(AW.E-E.AW), (17) 



8M 2 



^•(AW xE-Ex 



Vs = -8lP {A(2 '' <T ' B} ' (20) 
V " = -^4{ A<2, .°--( A(±)xE - ExA<±) )}. (21) 

1/7 = -i?s CT ' ExE - (22) 

and where the coefficients Cj are functions of aM and the running coupling a s in 
general. At tree level, their values are all unity. Note that the three-vector A^' 
refers to the spatial components of the covariant four- vector AR while E and B 
refer to the components E k and B k , respectively. 

The four-fermion contact interactions involving a quark and an antiquark x( x ) 
are given by: 

a 4 

^contact = -g 4 JpJ2{dl^( X )x(x)X f ( X )' l l J ( X ) + ( x ) CT%(x) • xK x )^{ x ) 
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+ ^(^(x)r X (x)x t (x)T>(x) +d 4 ^(x)T a ax(x) ■ X ] {x)T a a^{x))} . (23) 



a=l 



Analogous contact terms between pairs of quarks and pairs of antiquarks also occur 
but only affect baryons. 

The standard Euclidean gluon action is given by 

S G = a 4 J2£ G (x), (24) 

X 

where 

£g(x) = -4t E Tr [2 " JW*) - ^(*)] (25) 

is the usual single-plaquette lattice gauge field Lagrangian density. An improved 
gluonic action may be written 

S G = a^{\c G {x)- l -C^ 2 \x)\ (26) 

where Cq X2 ^ is comprised of 2 x 2 plaquette operators and is given by: 

4 2X2) (-) = 3^ E Tr[2 - fig(x) - fl« 2 > (*)], (27) 



with 



= 7 E ^a(^)^ r a(^+ae a )f/ /3 (a;+2ae a )f/ /3 (x + 2ae a + ae / 3) 

xL^_ a (x+2ae a +2ae / 3)L^_ a (x+ae a +2ae / g)C/_ / a(x+2ae ( g)C/_ / j(a; +aep), (28) 

and {(a, = {(n,v),(v,-n),(-n,-v),(-v,n)} for n ^ v. Presently, light 
quarks are neglected. 

The lattice NRQCD action is formulated in terms of the link variables U^(x) in 
order to preserve local gauge invariance, and as the lattice spacing a becomes small, 
it must tend to the action of the continuum theory. This can be shown at tree level 
in perturbation theory using the following relationship between the link variables and 
the gluon field G^(x) of the continuum theory: 

U^x) = exp iagA^x+^ip) ~ 1 + iagG^x), (29) 

where the lattice gluon field A^x) = G^x) + 0(a 2 ). Beyond tree level, however, one 
observes that large renormalizations are necessary to match the small a limit of the 
lattice action to the continuum form. These large renormalizations stem mainly from 
the higher-order powers of agA^ which occur in the expansion of U^. Such terms 
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generate ultraviolet divergences proportional to powers of a -1 and so are suppressed 
only by powers of g 2 and not a. 

A simple gauge-invariant procedure for improving the lattice NRQCD action by 
reducing the magnitudes of the renormalizations needed to reproduce the continuum 
theory has been suggested: replace all link variables that appear in the lattice 
action by U^/uq, where Uq is a parameter representing the mean value of the link [0. 
A gauge-invariant definition of this mean-field parameter may be written in terms of 
the mean plaquette: 

u = (^Tr^) 1 / 4 . (30) 

The parameter uq may be calculated using perturbation theory or may be measured in 
a simulation in order to include nonperturbative effects. The mean-field corrections 
introduced by this procedure are sometimes not small; for example, the cloverleaf 
magnetic and electric fields may be nearly doubled. Omission of such factors can 
result in significant underestimates of quantities, such as spin splittings, which depend 
on E and B. 

In lattice NRQCD simulations, a value for the lattice spacing a must be chosen 
subject to a few constraints. First, the cutoff A ~ n/a must be larger than the highest 
physically-relevant mass scale; thus, A ^> Mv is required. For bottomonium, v 2 ~ 0.1 
and M b ~ 5 GeV so that a < 2.0 GeV -1 ; for charmonium, v\ ~ 0.3 and M c ~ 1.7 GeV 
yielding a <C 3.4 GeV -1 . It is also desirable to choose A large enough so that 0(g 2 (A)) 
corrections to the coupling coefficients Cj and dj are small. However, NRQCD is a 
nonrenormalizable effective field theory and its coupling coefficients generally contain 
powers of A/M, causing problems if A is made too large. Furthermore, an overly 
large cutoff defeats the very purpose of NRQCD, which is to remove the scale M 
from the dynamics. Hence, it is important that A does not get much larger than M. 
Note that the requirement ai>r, where L is the extent of the lattice in number of 
sites and r is the root-mean-square radius of the system under study, can be easily 
satisfied in present-day simulations since ~ 2.4 — 5.0 GeV -1 for the ip family and 
r T ~ 1.0 — 3.6 GeV -1 for bottomonium. 

Once a suitable value for a is chosen, the bare quark mass M and the bare lattice 
coupling g must be fixed by reference to experiment. One very useful quantity for 
this purpose is the spin-averaged s—p splitting in quarkonium; due to an interplay 
between the Coulombic and linear forces, this splitting is essentially independent of 
aM. The masses of the lowest-lying ip and T mesons are also appropriate references. 
Of course, g and M run with a in such a way that low-energy physical predictions 
do not depend on the cutoff up to some order in v and a. In practice, it is much 
more convenient to first choose an appropriate value for g, and then fix a and M by 
comparison to experiment. 



7 



3 The Feynman Rules 



In a lattice gauge theory, the space of gauge transformations is finite so that gauge- 
fixing is not necessary. However, weak-coupling perturbation theory can only be 
applied if one fixes the gauge and extends the integration range of A b (x) using the 
familiar Faddeev- Popov technique ||. Hence, a gauge-fixing term S GF must be added 
to the NRQCD action. The Faddeev- Popov ghost action can then be determined from 
S GF 0. 

To facilitate the perturbative evaluation of scattering amplitudes, the gauge- 
invariant path integral measure must be expressed in terms of the lattice gluon field 

0: 



\xbfj, 



X 



where 



Sr, 



[VU\ = exp (S n 

2(1 - cos agA^x + |e M )) 



— > Tr In 
2 ^ 



(agA^{x + §e M )) 2 



(31) 



(32) 



and A b u c = if bcd A d u . To order g 2 , 



aV 
8 



E 

xb/i 



n 2 



(33) 



The Feynman rules are determined by expanding the total action in terms of the 
coupling g and Fourier transforming into momentum space. The coupling coefficients 
are written Cj = 1 + g 2 cf^ + 0(g 4 ). Rewriting the heavy quark-gluon action as 



S ( q ] = a 8 ^(y)G(y,x)^(x), 



xy 



then 



S, 



(n) 



d A k d A k' 



^(k')G(k';k)ilj{k), 



(2vr) 4 (2tt) 4 

where the Fourier transforms are defined on an infinite lattice by: 



ip(x) 
G(k';k) 



d A k 



-ik-x 



kn\<n/a (271-)' 

E 



jik' -y—ik-x 



G(y;x). 



(34) 
(35) 

(36) 
(37) 



xy 



The perturbative expansion of G(k'; k) in terms of the gluon fields takes the general 
form: 

G(k'; k) = (27^) 4 5 (4) (A; , - k) a" 1 ^*', k) 
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+ d(k' J k;q 1 ) A b u \( qi ) ^{k^k-q^M) 

+ E / d(k',k; qi ,q 2 ) Al\(q 1 )A%(q 2 )a^\k',t,q 1 ,u 1 ,b 1 -,q 2 ,u 2 ,b 2 



(3f 



+ 



where 



',k;qt,. . . ,q r ) 



n 



and ^4 M (gO is the Fourier transform of A^(x) defined by 

d 4 q 



(39) 



i=l 



(27T)" 



(40) 



Thus, the Feynman rules follow easily from the (dimensionless) functions. 



One of the simplest ways to compute the functions £^ is to first calculate the 



Fourier transforms of the basic operators which comprise G{k'\ k) and then combine 
these Fourier transforms appropriately. The perturbative expansion of the transform 
of each such operator will take the same general form as that for G(k'; k) given above. 
Of course, this general form also applies to transform products. Let 



C{ a py(k'; k) 



d 4 p 



A {a} (k';p) B m (p; k), 



(41) 



then the product rule for the £ functions may be expressed as follows: 

&^(V,b{<n,viMU) (42) 

r s r 

= 12^A {a} k ' + J2 ft! {&> Vi, bi} s i= ^i%~^(k - Qm, k] {qj, uj, bj} r j=s+1 ) 



1=1 



m=s+l 



Using this rule and the additivity of the £ functions, one can quickly build up the £^ 
functions. The Fourier transforms of all the necessary basic operators are presented 
in Appendix A. 

The lowest-order heavy quark propagator, shown in Fig. |l](a), is given by: 



Q%{k)=a^5 aP Q(k), 
where i, j are color indices and a,/3 are spin indices, and 



(43) 



Q(k) 



1 -e 



14 



K 2 (k) 2 \ 2 / « 2 (k) « 4 (k) K 2 (k) 2 \ 2n 



M 3 a 3 



nMa 3nMa f 2n 2 M 2 a 2 



(44) 



with K n (k) = J2 3 j= i sin n (kja/2). 
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In the Feynman gauge, the lowest-order gluon propagator, shown in Fig. ffl(b), is 
given by: 

(45) 



V^(k;\) = a 2 5 bc 5 riu D„(k;\) 



where 

D v (k;X) 



T ( E ) " ~, -s 2 (^) ( E ^(k a a) ) + a 2 A 



' fcjfyCL N 



a=l 



a=l 



(46) 



A small gluon mass A has been introduced to provide an infrared cutoff. If the simple 
gluon action in Eq. (|24j ) is used, the lowest-order gluon propagator is then 



V b T] c u (k;X) = a 2 5 bc 5 vu D(k;X), 



where 



D(k;X) 



<Esm 2 (^))+a 2 A 



-1 



(47) 



(48) 



The lowest-order vertex factors corresponding to interactions involving a heavy 
quark line and one to three gluons, shown in Fig. 0(c), may be written: 



Vx(k' 7 a, i; k, (3,j; q h »i, h) = -g (2tt) 4 ^ 4 ) (k'-k+qx) 



Y.<P T ij( i £ u (k',k-,q 1 ,K 



(49) 



V 2 (k',a,i;k,f],j;q 1 ,u 1 ,b l ;q2,U2,b2) = -ag 2 (27r)V 4) (k' - k + q 1 + q 2 ) 

^a^(r 6 n^) (^(fc'.^gn,^;?^^), (50) 

V 3 (k',a,i;k,f3,j;q 1 ,v 1 ,b 1 ;q 2 ,V2,b2;q3,V3,b3) = -a 2 # 3 (27r)V 4) (V- ^ + ^1 + ^2+^3) 

(t^t^t 6 ^) ..C0A;', fc; ?n,^ a ; ?^ 2 ; 

+^Tr(T^T^T^) C ^) (fc', fc; g n ,z/ Tl ; g T2 ,z/ T2 ; g T3 ,z/ r3 ) , (51) 



where P r is the group of permutations of r elements, a k are the standard Pauli spin 
matrices for k = 1,2,3 and a 4 is a 2 x 2 identity matrix in spin space. The (q q 
functions are obtained from the £q functions by neglecting the 0(g 2 ) corrections to 
the Cj coupling coefficients. These corrections show up in higher-order counterterms. 
Link variable renormalization in perturbation theory leads to the addition of the 
following order g 2 counterterm, shown in Fig. |l](d): 



V M0 (A;',CM;fc,/3,j) = ~4 2) -(2rr) 4 5^ (k'~ k)5 a p8, 
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/ « 2 (k) « 4 (k) K 2 (ky y- 1 / /. 2 (k)n 

V nMa 3nMa 2n 2 M 2 a 2 / \ M 3 a 3 y 

J A Mk) 5/t 2 (k) 2 \ / « 2 (k) «4(k) K 2 (k) 2 \ 

\ V M 3 a 3 M 3 a 3 / V nMa 3nMa 2n 2 M 2 a 2 j 

( ^(k) 2 \ /J7 4/c 2 (k) 4K 4 (k) 3« 2 (k) 2/t 2 (k) 2 \| 

V M 3 a 3 / \2Ma 3Ma 3Ma nM 2 a 2 nM 2 a 2 )y 



writing u = 1 + u g 2 + 0(g 4 ). 

The goal here is to determine the numerical values of the coupling coefficients Cj 
and dj and various renormalization factors for given values of the input parameters 
needed in lattice simulations, namely, the bare lattice coupling g and the bare heavy 
quark mass aM. Since these couplings and renormalization factors essentially absorb 
the relativistic effects arising from highly-ultraviolet processes, one expects that they 
may be calculated to a good approximation using weak-coupling perturbation theory. 
Using the above Feynman rules, the development of perturbative expansions for these 
quantities in terms of g is straightforward. However, there is no compelling reason 
to use the bare lattice coupling for the expansion parameter. In fact, recent work || 
suggests that g is a very poor choice of expansion parameter and that much better 
perturbation series result if one re-expresses the series in terms of a renormalized 
coupling g r defined in terms of some physical quantity and which runs with the 
relevant length scale. This is standard practice in continuum perturbation theory. Of 
course, if calculations could be carried out to all orders, then the choice of expansion 
parameter would be immaterial. 

To define a renormalized expansion parameter, a definition of the running coupling 
g r (fj,) and a procedure for determining the relevant mass scale fi must be given. A 
renormalization scheme || which defines the coupling such that the heavy quark 
potential has no gf or higher order corrections is particularly attractive. This scheme 
is physically motivated and produces O(g^) perturbative results in good agreement 
with simulation results for several different quantities. By absorbing the higher-order 
contributions to the heavy quark potential into g, 2 , it is hoped that higher-order 
contributions in other quantities will be small. The renormalized coupling g r ([i) 
is then given by the usual two-loop formula with A = 46.08A; at . The scale \i is 
determined by averaging lng 2 over the one-loop process of interest, where q is the 
loop momentum. Since the heavy quark parameters calculated here are ultraviolet 
divergent quantities, one expects fj, ~ iz/a. Here, g r shall be used to denote the value 
g r {fJ>) of the renormalized coupling at the appropriate scale. For example, at (3 — 5.7, 
g 2 (7i/a) ps 1.9 and at (3 — 6.0, g 2 {ir/a) ~ 1.7. Alternatively, g r could simply be 
added to the list of parameters which must be fixed in any simulation by reference to 
experiment. 
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4 The Heavy Quark Self-Energy 



The heavy quark self-energy may be defined by writing the inverse quark prop- 
agator Q~ l {p) in the form 

aG-\p)% = Q-Hp)S ij 5 aP - oEg,(p), (53) 

where i, j are color indices and a, j3 are spin indices. At order g 2 and v 4 and neglecting 
link variable renormalization for the moment, this self-energy is given by: 



ES»(P) = Km^ME^Cp; A) + A)), (54) 

where 

oEW(p;A) = ^a 4 £ / -^yQiP ~ k)D u (k; A) £,[cg> - fc,p; fc,i/)] 2 . (55) 

corresponding to the diagram in Fig. 0(a), and 

aE< B >(p; A) = --/Y.J ^D v (k;X)C^J P ,p-,k,u;-k,u), (56) 

corresponding to the tadpole diagram in Fig. ||(b). Note that e M = (—1, —1, —1, 1). 
The following properties of the ( functions are used to obtain the above results: 

(Vjk', k; k - k', u) = e, (Vjk, k'] k' - k, u), (57) 
C^{p,mh,v,-k,v) =0, (j = 1,2,3). (58) 

The self-energy is invariant under spatial reflections pj — > —pj, for j = 1,2,3, and 
transforms into its complex conjugate under p 4 — > —pi. 

In order to investigate E^(p;A) and £( B )(p;A) in the neighborhood of p = 0, 
the integrals in Eqs. (|55|) and (|56|) must first be evaluated. The usual initial step in 
the determination of such integrals is to use the change of variables z = exp(±zfc 4 a) 
to transform the integral over into a contour integral along the \z\ — 1 unit cir- 
cle. Unfortunately, the complicated pole structure of the vertex factors near z = 
makes difficult the evaluation of this contour integral by the residue theory. Due to 
this fact, the simplest procedure, approximating the four-dimensional integral by an 
appropriate summation as described below, is preferred. 

Because the integrand in Eq. ([56]) is a periodic analytic function of the real vari- 
ables with period 2tt when A > 0, T^ B \p\ A) is numerically well approximated by 
the discrete sum 

a^ B \ P - A) « -~ Y, £ D ^ WeLfan k > "i ~ k > ( 59 ) 

6 iV k u=l 
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In this sum, ak^ = li^n^j N where the n M take all integer values satisfying —N/2 < 
< N/2 for integral N. The error resulting from this approximation diminishes 
exponentially fast as iV — > oo. However, the rate of decay of this error is directly 
proportional to the mass gap aX, creating difficulties when aX is small. Fortunately, 
the decay rate of this error can be dramatically increased by making the following 
change of variables ||: — > — asm(k^) with < a < 1. This transformation 
maintains periodicity and effectively increases the mass gap so that the approximation 

aE< fl >(p; A) « -~ £ E g(k)D v (s(k); A)Cg 4 (p,p; s(k), v- -s(k), u), (60) 

where s At (fc) = k^ — asin(k^) and g(k) = n^=i[l — a cos(fc M )], converges much more 
quickly as N is increased. The parameter a should be chosen so as to maximize the 
effective mass gap: a = sech(w), where u satisfies aX ~ u — tanh(w). 

The above procedure is sufficient for evaluating (p; A) as long as the gluon mass 
aX is not set too small. However, the pole in the quark propagator is problematical 
when evaluating E^ A ^(p; A) near p = 0. To circumvent this, the contour for the a/c 4 
integral, which runs along the real axis from — it to it except near the pole, can 
be continuously deformed into a contour consisting of three line segments passing 
through the points — 7r — ► — 7r — iaX/2 — > 7r — iaX/2 — > ir. This contour is chosen for 
the following two reasons: for p ~ 0, the distance of closest approach to any pole is a 
maximum; and the contributions from the segments of the contour running parallel 
to the imaginary axis cancel due to the periodicity of the integrand. E^ A )(p ~ 0; A) 
can then be accurately obtained using the following approximation: 

oE< A >(p~0;A)«~2; E g(k)Q(p-r(k)) 

6ly k fl,U=l 

xD u (r(ky,X)e,[C^Jp-r(k), P] r(k),u)} 2 , (61) 

where r^(k) = k^ — asin(fc M ) — iaX5^ ll /2 and ak^ = 27m^/N with the taking 
all integer values satisfying —N/2 < < N/2. Also, a = sech(y) and aX/2 w 
y — tanh(y). In practice, the approximations in Eqs. ( |60"D and (^1|) are applied using 
increasing values of N until sufficient convergence is observed; typically, = 20 is 
adequate. Note that for values of p satisfying p x = p y = p z , the number of terms 
which must be independently evaluated in these sums may be dramatically reduced 
by exploiting the invariance of the summands under interchange of any two spatial 
components of k. 

The evaluation of £< A )(p ~ 0; A) and £< B >(p ~ 0; A) using Eqs. flBT] ) and flBTjp is 
feasible only if the gluon mass A is not set too small. However, the limits of these 
functions as A — > are actually required, necessitating the use of an extrapolation 
procedure. The expected behavior of these quantities as A tends to zero may be 
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expressed as an asymptotic expansion of the form 



Z(P; A) ~ E (€ ] + € ] ln« 2 A 2 ) (a\) m , (62) 

where the coefficient is known from the continuum theory since the infrared di- 
vergence is insensitive to the ultraviolet regulator. If one neglects the b$ terms 
for m > 0, then polynomial extrapolation using Neville's algorithm may be applied 
to /o(p; A) = E(p; A) — In a 2 A 2 . The extrapolation can be improved by apply- 
ing Neville's algorithm to the function f\(p; X) = fo(p',X) — aXd^ a x)^{p', A) + 2b^ 
which has no aXlna 2 X 2 term. Polynomial extrapolation of f'2{p',X) = /i(p;A) + 
(a 2 A 2 /2)(9 2 aA )E(p; A) + b^\ which does not suffer from a 2 A 2 lna 2 A 2 and aAlna 2 A 2 ef- 
fects, is another method. Since the gluon mass appears only in the gluon propagator, 
the derivatives of E(p; A) with respect to aX can be exactly and efficiently taken. 
In practice, polynomial extrapolation of all three functions /o(p;A), /i(p;A), and 
/2G9; A) is done using between six and twelve values of the gluon mass lying in the 
range 0.15 < aX < 0.8, and agreement among the results is verified. In calculat- 
ing £( B )(jo;A — > 0), one finds that only even powers of aX occur in its asymptotic 
expansion. In this case, the accuracy of the extrapolation can be increased by explic- 
itly excluding the odd powers in the extrapolating polynomial. Neville's algorithm is 
then applied to the functions fo(p', A), fi(p; A) = fo(p; A) — a 2 A 2 <9( a 2 A 2)E(p; X) + b^ and 

h(p', A) = flip] X) + (a 4 A 4 /2)(9 2 a 2 A 2)£(jo; A) + b^ /2. Uncertainties in the extrapolated 
values are estimated by examining the spread of values in the Neville table and by 
comparing the results obtained using different sets of gluon mass values and using 
the different functions described above. 

The small v expansion of the zeroth-order heavy quark inverse propagator, keeping 
only those terms which are suppressed by no more than v 2 relative to the leading 
terms, is given by: 

x • . P 2 a 2 pja 2 ip 4 p 2 a 3 (p 2 ) 2 a 4 (p 2 ) 2 a 4 

Q W--^ + ^ + ^ + ^^-^w-^W + ---' (63) 

recalling that p^ is of order v 2 . The on-mass-shell quark then satisfies the following 
dispersion relation: 

+ < 64 > 

where a>o(p) denotes the value of the fourth component of the contravariant quark 
four-momentum p^ at the pole of the zeroth-order propagator. This dispersion rela- 
tion agrees exactly with the continuum form from full QCD to this order in v. 

In the MS renormalization scheme, the inverse propagator to order g 2 for full 
QCD in continuum Minkowski space has the form 

(1 - C cont g 2 ) (V - M(l + 5 cont g 2 )) (65) 
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where 



C cont = I ^(-4+ln(M 2 //x 2 ) + 21n(M 2 /A 2 )), (66) 

Scant = -^(4 + 31n(/i 2 /M 2 )). (67) 
127T v J 

Note that \x is the mass scale introduced by dimensional regularization and A is the 
gluon mass regulating the infrared divergence. Thus, the sole effect of the order g 2 
corrections is to renormalize the quark field and mass. If lattice NRQCD is to repro- 
duce the low-energy physical predictions of full QCD, then the order g 2 corrections to 
the heavy quark propagator in lattice NRQCD must also do no more than renormalize 
the heavy quark field and mass to the appropriate order in v. 

Explicit calculation shows that the heavy quark self-energy in lattice NRQCD has 
a small v representation of the following form: 

2 2 

a£(p) « g 2 {n - ipiatlx + |^-0 2 + ...}, (68) 

retaining only radiative corrections to the lowest-order terms in v as specified by the 
power counting rules of Ref. M. The on- mass-shell quark now satisfies a dispersion 
relation given by 

wb(p) « i(-g 2 Q /a + ^ (l + g 2 n, - g 2 Q 2 ) + ...). (69) 

Defining M r = Z m M, where Z m = 1 — g 2 Qi+ g 2 Q 2 is the mass renormalization factor, 
and p4 = £>4 — ig 2 Qo/a, the inverse propagator for small v may be written: 

n-if\~r ( ■- , P 2 Q 2 , PW , ip4P 2 a 3 (p 2 ) 2 Q 4 (p 2 ) 2 Q 4 \ / 7n s 

where = 1 — g 2 (flo + ^i) is the wavefunction renormalization parameter. Thus, 
the addition of a counterterm which shifts the energy by an overall amount g 2 VLo/a 
is needed in order to match the low-energy physical predictions of lattice NRQCD 
with those of QCD. Alternatively, one could simply shift the energies obtained in 
simulations using the action in Eq. [14] by an amount g 2 Qo/a for each heavy quark. 

A more convenient set of renormalization parameters may be obtained by defining 
Z^ — 1 — g 2 C, Z m = 1 + g 2 B, and p 4 = p 4 + ig 2 A/a. The parameters A, B, and 
C can then be calculated using A = —Qo, B = Q 2 — ^i, and C = Qq + flj, where 
Vl = a£(0)/# 2 , fti = m<9 P4a £(0)/# 2 , and fi 2 = 2Ma 2 <9 p 2 a2 £(0)/# 2 . Due to the 
complexity of the NRQCD vertex factors, the derivatives in these expressions must 
be taken numerically. Four- and five-point formulas are applied in the differentiation; 
only points which satisfy p x = p y = p z are used since the self-energy can be computed 
much more quickly at such points. 
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Results for the energy shift parameter A and the mass renormalization parameter 
B are presented in Tables [1] and 0, respectively. Both A and B are gauge invariant 
and infrared finite. The wavefunction renormalization parameter C has an infrared 
divergence of the form — (lna 2 A 2 )/67r 2 as the gluon mass is taken to zero. This 
divergence is cancelled in physical quantities by an infrared divergence occurring in 
the quark-gluon vertex correction. Values for the infrared-finite portion of C are 
given in Table 0. In these tables, the contributions from the quark-gluon loop and 
tadpole diagrams are given separately. Results using the simple and the improved 
gluonic actions are presented. Values for the stability parameter n are chosen such 
that the pole in the quark propagator Q(p) falls on the same side of the real axis in 
the complex p 4 plane for all allowed values of p and tends to move farther away from 
this real axis as |p| increases to its allowed maximum. For aM > 2, n is set to unity; 
for 1 < aM < 2, n = 2 is used. These values of n ensure the stability of the evolution 
equation for the quark Green's function. 

Contributions from the quark-gluon loop graph to the energy shift parameter A 
are small and decrease in magnitude as aM decreases. At large aM, there is little 
difference between the shifts obtained using the simple and improved gluonic actions. 
The tadpole contributions are large and contain power-law divergences which grow 
in magnitude as aM is decreased. Since high-momentum modes are more strongly 
damped in the improved gluon propagator, the ultraviolet-divergent tadpole terms are 
appreciably smaller in the case of the improved gluon action. These total downward 
shifts in the energy are nearly the same as those obtained in Ref. which used a 
much simpler heavy quark action. Calculations of A using Cj = for j = 3, 4, 5, 6, 7 
reveal that contributions to this parameter from the spin-dependent interactions are 
small. 

Contributions to the mass renormalization parameter B from the quark-gluon 
loop diagram are very small and do not vary proportionately with aM. For all values 
of aM, there is little difference between the results obtained using the simple and 
the improved gluon actions. The tadpole contributions are again large, growing in 
magnitude as aM is decreased. The tadpole terms in the case of the improved action 
are slightly smaller and contributions to this parameter from the spin-dependent 
interactions are small. The total values for B obtained here are appreciably larger 
than the order g 2 corrections to the mass renormalization calculated in Ref. (see 
Ref. 0). 

The tadpole diagrams do not contribute to the heavy quark wavefunction renor- 
malization parameter C. The values for the infrared-finite portion of C are very 
small and become increasingly negative as aM is decreased. There is little difference 
between the results obtained using the simple and the improved gluonic actions. The 
magnitudes of the wavefunction renormalization corrections are much smaller than 
those obtained in Ref. |3j (see Ref. |10|1). 

The contributions to the energy shift and mass renormalization parameters from 
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the link variable renormalization counterterm are given by: 



SA 




7 



) 



(71) 



1 + 



2Ma 
1 



SB 



2 



3 



) 



(72) 



3 



AnMa 



+ 



M 2 a? 



No change in the wavefunction renormalization occurs. As previously stated, the pur- 
pose of link variable renormalization is to enhance the similarities between lattice and 
continuum gauge-field operators, especially those depending on the cloverleaf electric 
and magnetic fields. Consequently, the counterterms introduced by this renormaliza- 
tion offset the large tadpole contributions which commonly afflict lattice perturbation 
theory, offering a means of improving its convergence. As shown in Table |], the or- 
der g 2 corrections to the heavy quark renormalization parameters are small once the 

(2,) 

mean-field corrections are taken into account. In this table, the value u = —0.083 
obtained by evaluating Eq. ^ in perturbation theory for the simple gluonic action is 
used. For g 2 ~ 2 and a ~ lGeV -1 , the b quark receives approximately a 1% lowest- 
order correction to its mass and the c quark receives about a 7% mass correction. 

5 Conclusion 

Lattice NRQCD is an effective field theory which promises to make possible high- 
precision numerical studies of heavy quark systems. It is essentially a low-energy 
expansion of QCD in terms of the mean velocity of the heavy quarks in a typical 
heavy-quark hadron. To fully define lattice NRQCD, the coupling strengths of its 
interactions must be specified. These are determined by requiring that lattice NRQCD 
reproduces the low-energy physical results of continuum QCD. Since the role of these 
couplings is to absorb the relativistic effects arising from highly-ultraviolet QCD 
processes, one expects that they may be computed to a good approximation using 
perturbation theory, provided the quark mass M is large enough. 

The heavy quark self-energy in nonrelativistic lattice QCD was calculated to 0(a s ) 
in perturbation theory. An action which includes all spin-independent relativistic 
corrections to order v 2 and all spin-dependent corrections to order v 4 was used. The 
standard Wilson action and an improved multi-plaquette action were used for the 
gluons. Results for the mass and wavefunction renormalization and an overall energy 
shift were obtained. Contributions from the quark-gluon loop graph were found to 
be very small; however, the tadpole contributions were large. The values of these pa- 
rameters will be needed in future numerical simulations of quarkonium. The effective 
couplings will also be needed; calculation of these quantities is in progress. 

A tadpole improvement scheme in which all link variables are rescaled by a mean- 
field factor Uq was also applied in perturbation theory. The main purpose of this link 
variable renormalization was to enhance the similarities between lattice and contin- 
uum gauge-field operators, especially those depending on the cloverleaf electric and 
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magnetic fields. An important consequence of this scheme was a significant offsetting 
of the large tadpole contributions to the heavy quark renormalization parameters. Us- 
ing a perturbative approximation to the mean plaquette for uo, the tadpole- improved 
heavy quark renormalization parameters were shown to be small. This scheme offers 
a means of improving the convergence properties of lattice perturbation theory. 
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Appendix 

The momentum-space representations of various components of the lattice NRQCD 
action are presented in this Appendix. Below, ( functions are introduced and are de- 
fined by 

d r) 0'> k; gi,i/i,6il • • • I Qr^rA) = /C?^', k ; q^i, Wr) T bl ■ ■■T K . (73) 

Also note that 5^ is an "anti-Kronecker delta function" and is trivial when /x = v 
and unity otherwise. 



a uiiiiy umtsi wise. 

In momentum space, the link variable becomes 
U„(k'; k) = a 4 Y, e l(k '- k) - x e- tk » a exp 



iagA^x + |e M ) 



= (27r) 4 5 (4 \k> - A;) e — ( fc M+ fc M)/2 + mg J d ( jfe / > fc . g)e — K+^A^q) 
J d(k', k; Ql , q2)e- ia ^ +k ^ 2 AMA,(q 2 ) 
m "jl J rf(A;',A;;g 1 ,g 2 ,g3)e-^ +fc '' )/2 i M (g 1 )4(g 2 )i Ml 



a 2 g 2 



<? 3 ) + 



The momentum-space representation of may be written: 

CUk',k) = -ism[p' + kU 
C^wC^'^^i^i) = iacos[^(k' + k) fl ] 5^ U1 , 
( { *l(±)(k',k;q 1 ,v 1 ;q2,v 2 ) = i y sin[|(A;' + k)^ <Wi<W 2 , 

CI 

Cfl(±)(k',k;qi,u l ;q 2 ,u 2 ;q 3 ,iy 3 ) = -i— cos[-(fc' + fc)„] 5 ^5^5^ 



(74) 



(75) 
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The improved symmetric derivative has the following momentum-space representa- 
tion: 

C^( ±) (k', k) = i {-\ sin[^' + kU + \ sin[a(k' + *)„]) , (76) 
C^l(±)(k' ,k;Qi,^i) = *|<Wi (^4cos[^(A;' + A;) M ]-cos[a(A;' + A;) M ]cos(^gi^)^ , 
( { *l(±)(k',k;q 1 ,v 1 ;q 2 ,V2) = ia 2 5^ Ul 5^ U2 (J sin[|(fc' + 

-i sin[a(A; / + cos(^gi M ) cos(^ 2/1 ) - g cos^A;' + fc) M ] sin[|(gi - g 2 ) M ]) , 
C%±)(k',k; qi,vu q 2 ,v 2 ; qs^s) = ia 3 '5 ^5 ^ 2 S ^ 3 cos[^(k' + fc) M ] 

+ ^cos[a(k'+k)^\ cos[^(k' -k)^] + ^cos[a(k' + k + ^q 1 -^q 3 ) fl ) cos(^q 2fl ) 
The Fourier transform of the lattice Laplacian may be written: 

CU(k',k) = -^^(k'+ku 

Cl2A(2)(^'^;?i, i/ i) = 2asin[-(A;' + k) ui ] 5 4jl/1 , 
C^ > A (2){k',k;q l ,u l ;q 2 ,u 2 ) = -a 2 cos[^(k' + k) Ul ] 5 UuU2 5 A>Ul , (77) 

Co2aP) i k ' , k; q!,^; q 2 ,u 2 ; q 3 , u 3 ) = -— sm[-(k' + k) ui ] 5 UuU2 5 U2tU J 4tUl . 
Similarly for the improved Laplacian: 

Cw&> k) = -4 £ (sin 2 [^(A;' + fc),-] + ^sin 4 [^(A;' + fc),]) , (78) 
(^(^(^^^i^i) = 1*4,1^ ^8sin[|(fc / + fc^J-sin^^'+A;)^] cos(^gi^)) , 



Q2Aw( fc '' fc !9i^i; 92,^2) = a 5^,^*4,^ (-^cos[-(A; / + A;)^ 



+ i cos[a(/c / +A;) l/1 ] cos(^<7i„J cos^g^J- ^ sin[a(/c / + A;) i , 1 ] sin[|(gi-g 2 )^j , . 



(S > A(2 ) (k',k;q 1 ,is 1 ;q 2 ,is 2 ;q 3 ,is 3 ) = a?8 vul/2 8 V2 ^8 i>vl ( — - sin[-(fc' + k) Ul ] 



+ ^sm[a( y k' + k + ^q 1 -^q 3 ) Ul ] cos(^q 2u2 ) + ^ sm[a(k' + k) Ul ] cos[|(fc'- fc)ij) . 

Another important operator is the cloverleaf field strength tensor. Its Fourier 
transform is given by 

F, v (k';k) = a 4 Y,e l(k '~ k) - x F, u (x) 
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= W I d(k',k;q) [/» A b u (q) - f*(q) A^q)} T b 
b 

+igY: I d(k',k;q u q 2 ) T b ^T b 4 /£(<&,<&) A b J (q,) A b J (q 2 ) 



bib 2 



-/£(?*<fc)^(?i)^(^ 

+iag 2 fd(k',k; qi ,q 2 ,q 3 ) (t^T** - I^T^T 6 *) 



bl&2&3 



*»,<&) ^(gi)4 2 (?2)^ 3 (?3) + <&,<&) A b ^ qi )A b *{q 2 )A b *(qz) 

-/£(<Zi,52,<fe) A b J(q 1 )Al%q 2 )A b l ^(q 3 ) - f^( qi ,q 2 ,q 3 ) A* ( Ql ) A b * (q 2 ) A b J (q 3 ) 
-L E M,Q2,qi) A^( qi )A b ;(q 2 )A b ^q 3 ) + f^( qi ,q 2 ,q 3 ) A b J (qi)A b u 2 (q 2 )A b ^ (q 3 ) 

82, qi) A b ^ qi )A b *{q 2 )A b *{q 3 ) + f^( Ql , q 2 , q 3 ) i^i)^)^)} 

+0(g 3 ) (79) 



where 



fji(q) = sin(ag^)cos(^^), 



a 



2 



92) = cos[-(gi + g 2 ) M ] sin[-(gi + q 2 ) v ] sin[-( gi - g 2 ) M ], 



92) = ^[cos(^g lM ) cos(a9i„ + °^q 2v ) + cos(^g 2l ,) cos(ag 2/ , + 

+ cos{aq 2fl + ^q ljX ) cos{aq lu + ^q 2u ) - cos(^g^) cos(|g 2l /)], 

/^(9i,92,93) = cos[|(gi + g 2 + g3) M ]|^sin[a(g 1 + g 2 + g 3 ) ; ,] 

- cos [t|(<?i + ?2 + 93)1/] cos[^(gi - g 3 )„] sin(|g2i/)|, 

f^(Qi,Q2,q3) = cos[^(q 1 + q 2 + q 3 )^]cos[^(q 1 +q 2 + q 3 ) u }cos[^(g 1 + q 2 ) f ,} sm(^q 3 „) 

+ ^ cos[a(g 2 + ^g 3 )^] sin[|(gi + g 2 )i/], (80) 



92, 9s) = cos[|(gi + g 2 + g 3 ) M ] sin[a(^gi + g 2 + ig 3 ) J cos[^(gi - g 3 ) M ] 



Using = A b (-q), one can easily check that F b u (k'; k)* = F b u (-k'; -k). Fur- 

thermore, F^ v {k'] k) = —F u ^(k'; k), as required. Note the absence of a TrT bl T b2 term 
in the order g coefficient. 

The Fourier transform of the improved cloverleaf field strength tensor has the 
same form as that for F M „, but the functions f^ u , . . . must be replaced by the following 
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functions: 

= - cos ( a ftO - cos(aq u ))f^{q), 
/^(ft> ft) = ^{(5 - cos[a(gi + q 2 )^\ - cos[a(gi + g 2 )„])/* (gi, ft) 

-sin[a(gi + ^ft^Z/lXft) + sin[a(-gi + ft),,]./^^)}, 
/^(ft, ft) = ^{(5 - cos[a(g! + g 2 )J - cos[a(gi + ft)^])/^(ft, ft) 

+ sin[a(gi + -ft) i ,]/4(ft) + sin[a(-?i + ft) M ]/^(ft)}, 

/^(ft,ft,ft) = ^{(5-cos[a(gi+? 2 +? 3 ) /1 ]-cos[a(gi+? 2 +ft) I /])/2(ft,ft,ft) 

X 1 

+ sin[a(gi + ft + ^W^i, *0 - sin[a(-gi + q 2 + ft^l/^ft, ft) 

-i cos[a(g! + l -q 2 + ^ft) M ]/i(ft) - ^ cos[a(^gi + ^g 2 + ft)^]/i(ft) 
+ cos[a(^gi + q 2 + ^ftWi(ft)}, 

/J(ft>ft>ft) = ^{(5-cos[a(?i+g 2 +ft) /1 ]-cos[a(?i+g 2 +ft) i /])/^(ft,ft,ft) 

1 1 

- sin[a(gi + q 2 + -g 3 )/,]/^(gi, ft) - sin[a(-gi + q 2 + ft^J/^ft, ft) 

1 11 

-- cos[a(-gi + -q 2 + ft)z,]/^(ft)}, 

/^(ft,ft,ft) = ^{(5-cos[a(gi+? 2 +? 3 ) /1 ]- cos [ a (ft+ft+ft)^)/^(ft>ft>ft) 

- sin[a(gi + ft + ^3)^/^(51, ft) - sin[a(^gi + g 2 + ft) i ,]/^(ft, ft) 

- cos[a(ift + ft + ^ft),]/i(ft)}- (81) 
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Figures 

Figure 1: Various Feynman diagram elements. A curly line represents a gluon; a 
double solid line indicates a heavy quark, (a) Heavy quark propagator; (b) gluon 
propagator; (c) lowest-order vertices involving a heavy quark line and r gluons; (d) 
the 0(g 2 ) counterterm from link variable renormalization. 

Figure 2: Two Feynman diagrams which contribute to the heavy quark self-energy. 
A curly line denotes a gluon; a double solid line denotes a heavy quark. 
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Table 1: The energy shift parameter A for various values of the product of the bare 
heavy quark mass M and the lattice spacing a. The contribution to A from the quark- 
gluon loop diagram of Fig. [|(a) is denoted by Ai(A) for the improved gluon action of 
Eq. (|26| ) and by A S (A) for the simple gluon action of Eq. fl2"4|). The contribution from 
the tadpole diagram of Fig. ||(b) is denoted by Ai(B) and A S (B) for the improved 
and simple gluon actions, respectively. For aM > 2, the stability parameter n is set 
to unity; for 1 < aM < 2, n = 2 is used. Extrapolation uncertainties are no larger 
than ±0.0001. 
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Table 2: The heavy quark mass renormalization parameter B for various values of the 
product of the bare heavy quark mass M and the lattice spacing a. The contribution 
to B from the quark-gluon loop diagram of Fig. ^|(a) is denoted by Bi(A) for the 
improved gluon action of Eq. (|26"D and by B S (A) for the simple gluon action of Eq. fl24|). 
The contribution from the tadpole diagram of Fig. |](b) is denoted by B,i{B) and 
B S (B) for the improved and simple gluon actions, respectively. For aM > 2, the 
stability parameter n is set to unity; for 1 < aM < 2, n = 2 is used. Extrapolation 
uncertainties are no larger than ±0.0001 for the tadpole values and ±0.0002 for the 
contributions from the quark-gluon loop diagram. 
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Table 3: The infrared-finite portion of the heavy quark wavefunction renormalization 
parameter C for various values of the product of the bare heavy quark mass M and 
the lattice spacing a. The infrared-finite contribution to C from the quark-gluon loop 
diagram of Fig. |2](a) is denoted by Cj for the improved gluon action of Eq. ( p6|) and 
by C s for the simple gluon action of Eq. (|24D . The tadpole diagram of Fig. ||(b) does 
not contribute to this parameter. For aM > 2, the stability parameter n is set to 
unity; for 1 < aM < 2, n = 2 is used. Extrapolation uncertainties are no larger than 
±0.0001. 
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Table 4: Tadpole improvement of the energy shift parameter A and mass renormaliza- 
tion parameter B for various values of the product of the bare mass M and the lattice 
spacing a. Results are given for the simple gluon action only. The renormalization 
parameters without tadpole improvement are denoted by A s and B s ; the improved 
parameters are denoted by A s and B s . The perturbative value iiif 1 = —0.083 is used 
for the mean-field parameter. 
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5.00 
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